Sampling methods, reconstruction methods and devices for sampling and/or reconstructing multidimensional signals

ABSTRACT

Reconstruction method and devices for two-dimensional signals that are not bandlimited but have a parametric representation with a finite number of degrees of freedom. The signal is reconstructed from the samples obtained after a suitable filtering with a smoothing kernel. Projection techniques allow reducing the problem to a combination of one-dimensional subproblems.

[0001] Sampling methods, reconstruction methods and devices for samplingand/or reconstructing multidimensional signals Sampling theory hasexperienced a strong research revival over the past decade, which led torefinement of original Shannon's theory and development of more advancedformulations with direct impact in signal processing and communications.For example, the traditional sampling paradigm for representation ofbandlimited functions can be extended to classes of signals, such assplines, that are not bandlimited but live on a subspace spanned by agenerating function and its shifts.

[0002] International Patent Application WO 02/078197, in the name of theapplicant and which is hereby incorporated by reference, developssampling schemes for a larger class of non-bandlimited signals, such asstream of Diracs, non-uniform splines and piecewise polynomials. Acommon feature of these signals is that they have a parametricrepresentation with a finite number of degrees of freedom, and can beperfectly reconstructed from a finite set of samples.

[0003] On the other hand, while there are many interesting results forone-dimensional signals, the problem becomes more involved when going tohigher dimensions. Furthermore, most of the multidimensional samplingalgorithms encountered in practice still rely on results from abandlimited case, which may lead to unnecessarily high computationalload, particularly for those classes of signals that could presumably berepresented by a finite number of samples. Therefore, a very interestingbut challenging question is whether it is possible to come up withpractical methods for sampling multidimensional signals having finitecomplexity, which would allow for perfect reconstruction from a finiteset of samples.

[0004] One aim of the present invention is to provide sampling andreconstruction methods for multidimensional signals that overcome thelimitations of the prior art.

[0005] Another aim of the present invention is to provide sampling andreconstruction methods for multidimensional signals that decompose theproblem into a number of one-dimensional subproblems.

[0006] Another aim of the present invention is to provide sampling andreconstruction methods and devices for determining exactly the directionand the position of details in a set of images, or for determiningexactly the position of photo cameras from the images taken therefrom.

[0007] The above aims are attained by the object of the appended claims.

[0008] Various embodiment of the sampling and reconstruction method ofthe invention applied to different classes of bidimensional signals withfinite rate of innovation are explained in the specification. The oneskilled in the art will understand however that the invention is notlimited to the specific examples given, and that the advantageoustechnical effects can be obtained by sampling and reconstructing otherkinds of signals with a finite rate of innovation.

[0009] The invention will be better understood with the help of thefigures in which:

[0010]FIG. 1 diagrammatically illustrates a device and method ofsmoothing and sampling a signal;

[0011]FIG. 2 represents the projection of some represents the projectionof the set of three Diracs onto one line according to one aspect of theinvention

[0012]FIGS. 3a and 3 b represents the set of projecting lines along fourdifferent directions, obtained from the samples. Points where exactlyfour projecting lines intersect correspond to the Diracs in the set;

[0013]FIGS. 4a and 4 b refer to an embodiment of the inventionconcerning the reconstruction of the position of a point from a set ofimages

[0014]FIG. 5 refers to a further embodiment of the invention;

[0015]FIG. 6a represents a signal made up of 17 weighted Diracs;

[0016]FIG. 6b represents a 2-D sinc smoothing kernel;

[0017]FIG. 6c represents an approximation obtained by convolving thesignal of FIG. 6a with the kernel of FIG. 6b;

[0018]FIG. 6d represents the signal reconstructed according to oneaspect of the invention.

[0019] In a first embodiment of the present invention, we exploit theproperties of, for example, the Radon transform of multi-dimensionalsignals, and show that by taking a finite number of “filtered” lineintegrals, the problem can be reduced to its one-dimensional equivalent,which is much more convenient for algorithmic implementation.

[0020] The International Patent Application WO 02/078197 defines theconcept of rate of innovation of a signal, and describes sampling andreconstruction schemes for signals that are not bandlimited, but exhibita finite rate of innovation or a locally finite rate of innovation.

[0021] The rate of innovation of a one-dimensional signal is defined asthe number of degrees of freedom of the signal within a time period. Forexample a periodic Poisson process signal made of K weighted Diracpulses δ(t) and having period τ, can be fully specified by the Koccurrence times or shifts t_(k) and the K amplitudes or weights c_(k)of the Diracs. The signal can be written as $\begin{matrix}{{x(t)} = {\sum\limits_{n}^{\quad}\quad {\sum\limits_{k = 0}^{K - 1}\quad {c_{k}{\delta \left( {t - \left( {t_{k} + {n\quad \tau}} \right)} \right)}}}}} & (1)\end{matrix}$

[0022] wherein the degree of freedom of this signal is obviously 2K. Itsrate of innovation is thus finite and equal to 2K/τ.

[0023] On the other hand the signal (1) is not band-limited. Therefore,according to the classic sampling theorem of Whittaker, Kotelnikov andShannon, it can not be faithfully reconstructed by a finite number ofits samples. The International Patent Application WO 02/078197 teacheshow to reconstruct exactly non-bandlimited signals from a finite numberof samples y[n] of a smoothed signal y(t), obtained by convolution witha suitable smoothing kernel φ. The process is shown on FIG. 1.

[0024] The main purpose for introducing the rate of innovation p, isthat it can be directly related to the minimum sampling rate, or theminimum number of samples that can lead to perfect reconstruction.

[0025] The Continuous-time Fourier series (CTFS) coefficients of x(t)are defined by $\begin{matrix}{{X\lbrack m\rbrack} = {\frac{1}{\tau}{\sum\limits^{\quad}\quad {c_{k}^{{- }\quad 2\pi \quad t_{k}{m/\tau}}}}}} & (2)\end{matrix}$

[0026] If the signal x(t) is convolved with a sinc smoothing kernel ofbandwidth [−K/τ, K/τ] then we obtain the lowpass approximation y(t)given by $\begin{matrix}{{y(t)} = {\sum\limits_{m = {- K}}^{K}\quad {{X\lbrack m\rbrack}^{\quad 2\pi \quad t\quad {m/t}}}}} & (3)\end{matrix}$

[0027] Suppose that the lowpass approximation y(t) is sampled atmultiples of T, we obtain τ/T discrete samples defined by$\begin{matrix}{{{y_{s}\lbrack l\rbrack} = {{y({lT})} = {{\sum\limits_{m = {- K}}^{K}\quad {{X\lbrack m\rbrack}^{\quad 2\quad \pi \quad m\quad {{lT}/\tau}}\quad l}} = 0}}},{{\ldots \quad {\tau/T}} - 1}} & (4)\end{matrix}$

[0028] as long as the number of samples is larger than the number ofvalues in the spectral support of the lowpass signal y(t), that isτ/T≧2K+1, the values of the 2K+1 central values of X[m] coefficients isdetermined by relation (4). The X[m] can be obtained by applying theFourier transform to the samples y_(s)[l] of the filtered signal, in awell known way.

[0029] The locations and weights of all the Diracs in (2) can becompletely determined by 2K+1 CTFS coefficients X[m]. From (2) we havethat the CTFS coefficients X[m] are linear combinations of K complexexponentials. $\begin{matrix}{u_{k}^{m} = {{\exp \left( {{- 1}\frac{2\pi \quad {mt}_{k}}{\tau}} \right)}\quad {with}\quad {weights}\quad \frac{1}{\tau}c_{k}}} & (5)\end{matrix}$

[0030] Thus to find the locations tk we may for example proceed byfinding the annihilating filter H(z) with coefficients 1, H[1], H[2], .. . H[K] defined as

H(z)=1+H[1]z ⁻¹ +H[2]z ⁻² + . . . +H[k]z ^(−k)  (6)

[0031] which satisfies $\begin{matrix}{{{\sum\limits_{k = 0}^{K}\quad {{H\lbrack k\rbrack}{X\left\lbrack {m - k} \right\rbrack}}} = {{0\quad m} = 1}},{\ldots \quad K}} & (7)\end{matrix}$

[0032] It can be shown that the k shifts t_(k) of the Dirac are givenby, or at least can be retrieved from the zeros of H(z).

[0033] The filter H(z) can be found by solving the following structuredlinear equation system: $\begin{matrix}{{\begin{bmatrix}{X\lbrack 0\rbrack} & {X\left\lbrack {- 1} \right\rbrack} & \cdots & {X\left\lbrack {{- K} - 1} \right\rbrack} \\{X\lbrack 1\rbrack} & {X\lbrack 0\rbrack} & \quad & {X\left\lbrack {{- K} - 2} \right\rbrack} \\\vdots & \quad & \quad & \vdots \\{X\left\lbrack {K - 1} \right\rbrack} & {X\left\lbrack {K - 2} \right\rbrack} & \cdots & {X\lbrack 0\rbrack}\end{bmatrix} \cdot \begin{bmatrix}{H\lbrack 1\rbrack} \\{H\lbrack 2\rbrack} \\\vdots \\{H\lbrack K\rbrack}\end{bmatrix}} = {- \begin{bmatrix}{X\lbrack 1\rbrack} \\{X\lbrack 2\rbrack} \\\vdots \\{X\lbrack K\rbrack}\end{bmatrix}}} & (8)\end{matrix}$

[0034] Because the matrix is a Toeplitz system, fast algorithms areavailable for finding the solution. This system has an unique solutionif the matrix is invertible, which is the case if all K Diracs in x(t)are distinct.

[0035] Given the coefficients H, the filter H(z) can be factored intoits roots $\begin{matrix}{{H(z)} = {\prod\limits_{k = 0}^{K - 1}\quad \left( {1 - {u_{k}z^{- 1}}} \right)}} & (9)\end{matrix}$

[0036] which leads to the shifts t_(k) using the above relation betweenu_(k) and t_(k).

[0037] Given the shifts t_(k) the K values can be found by solving$\begin{matrix}{{X\lbrack m\rbrack} = {\frac{1}{\tau}{\sum\limits_{k = 0}^{K - 1}\quad {c_{k}^{{- }\frac{\pi \quad m\quad t_{k}}{\tau}}}}}} & (10)\end{matrix}$

[0038] which leads to the following Vandermonde system $\begin{matrix}{\begin{bmatrix}{X\lbrack 0\rbrack} \\{X\lbrack 1\rbrack} \\\vdots \\{X\left\lbrack {K - 1} \right\rbrack}\end{bmatrix} = {{\frac{1}{\tau}\begin{bmatrix}1 & 1 & \cdots & 1 \\u_{0} & u_{1} & \cdots & u_{K - 1} \\\quad & \quad & \quad & \vdots \\u_{0}^{K - 1} & u_{1}^{K - 1} & \cdots & u_{K - 1}^{K - 1}\end{bmatrix}} \cdot \begin{bmatrix}c_{0} \\c_{1} \\\vdots \\c_{K - 1}\end{bmatrix}}} & (11)\end{matrix}$

[0039] which can be solved in a known way.

[0040] Similar reconstruction methods are possible for other classes offunctions having a finite rate of innovation, namely bilevel signals,splines, polynomial signals, piecewise linear signals, piecewisepolynomial signals and approximations thereof. Moreover the smoothingkernel can be other than a sinc function, for example a Gaussianfunction, a spline function, a box function or a hat function, as it isknown from PCT/EP02/03380.

[0041] The rate of innovation of a multi-dimensional signal can bedefined in a perfectly analogous fashion. Consider for example atwo-dimensional bandlimited signal g(x,y), with the Fourier transformthat is nonzero over a finite region R in the frequency space(f_(x),f_(y)). If we let 2B_(x) and 2B_(y) represent the widths in thef_(x) and f_(y) directions of the smallest rectangle that encloses theregion R, According to the sampling theorem of Whittaker, Kotelnikov andShannon, the signal can be perfectly represented by the uniformly spacedsamples c_(min) via convolution with the sinc kernel. $\begin{matrix}{{g\left( {x,y} \right)} = {\sum\limits_{m = {- \infty}}^{+ \infty}{\sum\limits_{n = {- \infty}}^{+ \infty}{c_{m\quad n}\frac{{\sin \left( {{\pi \left( {x - {m\quad X}} \right)}/X} \right)} \cdot {\sin \left( {{\pi \left( {y - {n\quad Y}} \right)}/Y} \right)}}{{{\pi \left( {x - {m\quad X}} \right)}/X} \cdot {{\pi \left( {y - {n\quad Y}} \right)}/Y}}}}}} & (12)\end{matrix}$

[0042] where X and Y are such that X≦1/2B_(x) and Y≦1/2B_(y), andg_(min)=g(mX, nY).

[0043] The above result can be interpreted in the following way. Thebandlimited signal g(x,y) can be considered as having 1/X and 1/Ydegrees of freedom per unit of length in x andy directions respectively,which we call the rates of innovation ρ_(x) and ρ_(y). The total rate ofinnovation of the signal is then defined as the sum of the rates ofinnovation for each independent variable ρ=ρ_(x)+ρ_(y).

[0044] A more general form of the above expression for g(x,y), whichextends the subspace of bandlimited functions, is obtained by replacingthe function sinc(x,y) by a more general template, the so-calledgenerating function φ(x,y). $\begin{matrix}{{g\left( {x,y} \right)} = {\sum\limits_{m = {- \infty}}^{+ \infty}{\sum\limits_{n = {- \infty}}^{+ \infty}{c_{m\quad n}{\phi \left( {\frac{x - x_{m}}{X},\frac{y - y_{n}}{Y}} \right)}}}}} & (13)\end{matrix}$

[0045] x_(m) and y_(n) being arbitrary shifts. Clearly this substitutiondoes not alter the rate of innovation ρ of g(x,y). When φ(x,y)=δ(x,y)and both x_(n)−x_(n−1) and y_(n)−y_(n−1) are i.i.d random variables withexponential density, then g(x,y) describes a 2-D Poisson process.Examples of other classes include simple lines and polygonal lines,planar parametric curves, as well as bi-level signals whose boundarieshave a finite number of degrees of freedom.

[0046] We will continue the analysis by considering a signal g(x,y) thatis made up of M Diracs generated by a 2-D Poisson process, which has asimple parametric representation, but whose algebraic structure providesa good insight into the fundamental principles inherent in algorithmsfor more complicated signals. The concept we present is based onsampling the Radon transform of the signal g(x,y), and offers apossibility of decomposing the problem into a set of 1-D equivalents.

[0047] Let the signal g(x,y) be represented as $\begin{matrix}{{g\left( {x,y} \right)} = {\sum\limits_{k = 0}^{K - 1}{c_{k}{\delta \left( {{x - x_{k}},{y - y_{k}}} \right)}}}} & (14)\end{matrix}$

[0048] and let R_(g)(p,θ) be the Radon transform of g(x,y), defined by

R _(g)(p,θ)=∫g(x,y)δ(p−x·cos(θ)−y·sin(θ)) dxdy  (15)

[0049] On FIG. 2 a two-dimensional signal g(x,y) composed of three Diracpulses is represented on a xy plane. The Diracs are at the positions(x₁,y₁), (x₂,y₂), (x₃,y₃). For given angle θ₀ and distance p from theorigin, the Radon transform operation R_(g)(p,θ₀) is given by theintegral of g(x,y) along the line l_(p,θ) defined byp=x·cos(θ)+y·sin(θ), crossing the axis x at an angle θ₀ and distant pfrom the origin O.

[0050] For any given angle θ₀, the Radon transform R_(g)(p, θ₀) can beinterpreted graphically as the one-dimensional projection of thefunction g(x,y) onto the line l′_(θ) passing from the origin andorthogonal to the line l_(p,θ) a as it is visible on FIG. 2.

[0051] In the case above, the Radon transform R_(g)(p,θ₀) of the signalg(x,y) composed of three Dirac pulses will be also composed of threeDirac pulses, located at the points p₀₁, p₀₂, p₀₃, corresponding to theprojections of the 2-dimensional Diracs at (x₁,y₁), (x₂,y₂), (x₃,y₃),along the projection lines 71, 72, 73 respectively, on the line l′_(θ).

[0052] The Radon transform R_(g)(p,θ₀) possess therefore a finite rateof innovation, and can be represented as a one-dimensional weighted sumof M₀ Diracs $\begin{matrix}{{R_{g}\left( {p,\theta_{0}} \right)} = {\sum\limits_{k = 0}^{M_{0} - 1}\quad {a_{0\quad k}{\delta \left( {p - p_{0k}} \right)}}}} & (16)\end{matrix}$

[0053] In the general case, for almost every value of θ₀, each Diracwill project a distinct counterpart on the line l′_(θ), thus M₀=M almostalways. For some particular choice of θ₀, however, there will be morethan one Dirac spike on the path of integration, therefore in generalM₀≧M holds.

[0054] This fact is exploited in this embodiment of the invention totackle the problem in the 2-D case. Since the signal consisting of M 1-DDiracs can be perfectly recovered from a set of 2M samples, as it isknown from the above mentioned International Patent Application WO02/078197, by using either the sinc or the Gaussian sampling kernel, wecan take advantage of that result and develop a sampling scheme for 2-Dsignals as well. Namely, instead of taking the line integral in theequation (4), we can replace the δ-function with the appropriate kernel,such as the sinc function. In other words, we can consider the filteredprojection of the signal g(x,y) $\begin{matrix}{{{\overset{\_}{R}}_{g}\left( {p,\theta} \right)} = {\int{{g\left( {x,y} \right)}B_{p}\frac{\sin \left( {\pi \quad {B_{p}\left( {p - p_{\theta}} \right)}} \right)}{\pi \quad {B_{p}\left( {p - p_{\theta}} \right)}}{x}{y}}}} & (17)\end{matrix}$

[0055] where p_(δ)=x·cos(θ)+y·sin(θ). An interesting property emergingform this formulation is that for any angle θ₀, the projection of g(x,y)becomes a convolution of its Radon transform R_(g)(p,θ₀), whichrepresents a 1-D stream of Diracs, and the sine sampling kernel, i.e.$\begin{matrix}{{{\overset{\_}{R}}_{g}\left( {p,\theta} \right)} = {{R_{g}\left( {p,\theta_{0}} \right)}*B_{p}\frac{\sin \left( {\pi \quad B_{p}p} \right)}{\pi \quad B_{p}p}}} & (18)\end{matrix}$

[0056] The above equation implies, by virtue of the method disclosed inInternational Patent Application WO 02/078197, that the locationspk andweights a_(k) of the 1-D Diracs, defined by (16), can be obtained fromN≧2M₀ samples of {overscore (R)}_(g)(p,θ), that is,

p _(n)(θ₀)={overscore (R)} _(g)((p−nT _(p),θ₀) n=0,1, . . . N−1  (19)

[0057] where T_(p)=1/B_(p).

[0058] While the set of locations p_(0k) does not itself define g(x,y),the projection of g(x,y) onto M+1 straight lines entirely specifies thesignal. Assume therefore that we do the projection of g(x,y) onto M+1lines with different slopes, determined by angles θ₀, θ₁, . . . , θ_(M),as shown in FIGS. 3a and 3 b. By using the method described above, wecan solve for the coordinates p_(mk) and weights a_(mk), m=0,1, . . . Mof the 1-D Dirac streams along each line, and thus uniquely specify theset of “projecting” lines that extend perpendicularly from each p_(mk),as it is visible on FIG. 3b. Clearly, for any point that belongs to theset of 2-D Diracs, exactly M+1 projecting lines must intersect. Areverse statement, that the points where M+1 projecting lines intersectmust belong to the set, can be proved by counterexample. Namely, assumethat exactly M+1 such lines intersect at some point A. If A doesn'tbelong to the set of Diracs, then there must be at least one point fromthe set located at each of the lines l_(p) _(—) _(mk,θ) _(—) _(m),m=0,1, . . . M, which implies that g(x,y) is being made up of at leastM+1 Diracs, and that obviously contradicts our basic assumption.

[0059] The weights c_(k) can be found by solving the system of M linearequations that we can choose out of (M+1)2M equations (available fromthe set of projections of g(x,y)). Thus the signal g(x,y) composed of afinite set of M weighted 2-D Diracs can be perfectly reconstructed from2M(M+1) samples of the filtered Radon transforms {overscore(R)}_(g)(p,θ) defined in equation (17).

[0060] As can be seen from the above derivation, in a general case thealgorithm yields a unique solution by taking on the order of M² samplesof the Radon transform. We will show that in some cases of morecomplicated signals, the algorithm's complexity remains either the sameor can be even reduced, depending on the specific signal structure, withthe only difference in the sampling scheme being the use of analternative sampling kernel.

[0061] It is to be noted that in a physical set up, very often the lineintegral present in the Radon transform cannot be obtained orimplemented exactly. In a typical application the signal g(x,y)represents an image containing sharp details or point sources.Conventional image-recording devices are not capable to capture theoriginal signal. Rather, a smoothed version is obtained.

[0062] In a real optical system, for example, the image will alwayssuffer degradation from optics' finite resolution or atmosphericeffects. Likewise the projections obtained in an X-ray scan will bedegraded because the X-rays can not be collimated in an infinitelynarrow pencil beam. That is, effectively, a filtered or smoothed Radontransform is obtained.

[0063] While mathematically, the filtering and Radon transform can beinterchanged, in the physical set up one can only obtain the filteredprojection rather than the ideal Radon transform. In the case of theoptical system cited above, one has no direct access to the originalsignal g(x,y), but only to the smoothed image provided by the opticalsystem. In the case of an X-ray scan, projection and smoothing arecoupled together and can not be performed independently.

[0064] The degradations induced by the physical setup can in general bemodeled by an appropriate choice of the smoothing kernel; for example bya Gaussian kernel.

[0065] The reconstruction method of the present invention provides aperfect reconstruction of the original signal g(x,y) from a finite setof samples of the filtered Radon transform.

[0066] In particular, if the original signal g(x,y) represents an image,the inventions allow reconstructing a “superresolution” image whosesharpness is much superior to that provided by conventionalimage-recording apparatuses. It must be noted, however that the presentinvention is not limited to image treatment, and rather it consists of ageneral method or reconstructing a multidimensional signal from a set ofsamples.

[0067] The Radon transform has been chosen in the above example merelyto illustrate a particular mode of realization of the invention. Theinvention however includes also those algorithms in which variations orextensions of the Radon transform, or other transformations, are used toproject the multidimensional function onto an ensemble ofone-dimensional functions.

[0068] The sinc kernel employed in the above method does not constitutea limitative feature of the invention, but merely an example of one outof the many kernels utilizable in the present invention, for examplegaussian, splines or hat kernels could advantageously be used whenfast-decaying tails or a compact support are desirable.

[0069] In a second embodiment of the present invention, we deal withanother class of non bandlimited 2-D functions of finite rate ofinnovation, namely that of bi-level polygons.

[0070] Consider a signal g(x,y) that is a bi-level polygon uniquelydefined by its vertices located at points (x_(j),y_(i)), i.j=1 . . . M.Clearly g(x,y) has 2M degrees of freedom. If we now project g(x,y) on astraight line, as in the previous embodiment, the result will be a 1-Dpiecewise linear signal, we can therefore use the same approach andincorporating the sampling scheme for such 1-D signal into ouralgorithm, by replacing the δfunction in the Radon transform into anappropriate smoothing kernel.

[0071] A 1-D piecewise linear signal having M pieces can be uniquelyrepresented by 2M samples <f(p), φ⁽²⁾ (p−nT_(p))> of its convolutionwith the second derivative sinc function, here defined as an inverseFourier transform $\begin{matrix}{{\phi^{(2)}(p)} = {{IFT}\left\{ {\left( {j\quad \omega} \right)^{2}{{rect}\left( \frac{\omega}{2\quad \pi \quad B_{p}} \right)}} \right)}} & (20)\end{matrix}$

[0072] with B_(p)=1/T_(p). Due to the associativity of the convolutionoperator, the convolution f(p)* φ⁽²⁾ is equivalent to the convolution ofthe second derivative of the signal with the sinc kernel, f⁽²⁾ (p)*φ,thus the problem can be reduced to the one we have already analyzed,because the second derivative of a piecewise linear signal is again astream of Dirac pulses.

[0073] In other words the sampling scheme in the 2-D case should bemodified such that the δ-function in (14) is replaced by φ⁽²⁾. In thatcase taking at most 2M samples from each of the M+1 projections andapplying the same techniques will uniquely and exactly specify thecoordinates of the vertices.

[0074] In a further embodiment of the invention, a scheme is devised forsampling and reconstructing a bi-level signal with a piecewisepolynomial boundary: $\begin{matrix}{{g\left( {x,y} \right)} = \left\{ \begin{matrix}1 & {x \leq {p(x)}} \\0 & {x > {p(x)}}\end{matrix} \right.} & (21)\end{matrix}$

[0075] where p(x) Is a piecewise polynomial having Mpieces of maximumdegree R. We can directly apply the result from the 1-D case, and takethe samples of the projection of g(x, y) on the x axis, by using the(R+1)th derivative sinc kernel. Since the (R+1)th derivativep^((R+1))(x) is a collection of at most M weighted Dirac pulses, we needto take only 2M samples in order to exactly recover the signal.

[0076] In a further embodiment of the invention, a method is devised toreconstruct a superposition of 2-D Diracs, without resorting to theprojection into a number of one-dimensional subproblems.

[0077] Let g(xy) be a 2-D periodic signal of period T along bothdirections $\begin{matrix}{{g\left( {x,y} \right)} = {\sum\limits_{p,q}{\sum\limits_{k = 0}^{M - 1}\quad {a_{k}{\delta \left( {{x - {p\quad T} - x_{k}},{y - {q\quad T} - y_{k}}} \right)}}}}} & (22)\end{matrix}$

[0078] The Fourier series representation of g(x,y) $\begin{matrix}{{g\left( {x,y} \right)} = {\sum\limits_{m = {- \infty}}^{+ \infty}{\sum\limits_{n = {- \infty}}^{+ \infty}{{G\left\lbrack {m,n} \right\rbrack}{\exp \left( {\quad m\quad \omega_{0}x} \right)}{\exp \left( {\quad n\quad \omega_{0}y} \right)}}}}} & (23)\end{matrix}$

[0079] where ω₀=2π/T and the Fourier coefficients G[m,n] are given, asit is well known, by $\begin{matrix}{{G\left\lbrack {m,n} \right\rbrack} = {\sum\limits_{k = 0}^{M - 1}{a_{k}\exp \quad \left( {{- }\quad m\quad \omega_{0}x} \right){\exp \left( {{- }\quad n\quad \omega_{0}y} \right)}}}} & (24)\end{matrix}$

[0080] that is a linear combination of M complex exponentials. Inanalogy with the one-dimensional case, the knowledge of a finite numberof Fourier coefficients G[m,n] uniquely determines the positions (x_(k),y_(k)) and weights a_(k) in (22).

[0081] This can be done operatively by applying known singular valuedecomposition methods and subspace methods.

[0082] To make the notation simpler (24) can be written as$\begin{matrix}{{G\left\lbrack {m,n} \right\rbrack} = {\sum\limits_{k = 0}^{M - 1}{a_{k}w_{k}^{m}z_{k}^{n}}}} & (25)\end{matrix}$

[0083] ir, in matrix notation G=WAZ with W, A and Z defined as$\begin{matrix}{W = \begin{bmatrix}1 & 1 & \cdots & 1 \\w_{1} & w_{2} & \quad & w_{M} \\\vdots & \quad & \quad & \quad \\w_{1}^{P - 1} & w_{2}^{P - 1} & \quad & w_{M}^{P - 1}\end{bmatrix}} & (26)\end{matrix}$

 A=diag(a ₁ ,a ₂ , . . . ,a _(M))  (27) $\begin{matrix}{Z = \begin{bmatrix}1 & z_{1} & \cdots & z_{1}^{Q - 1} \\1 & z_{2} & \quad & z_{2}^{Q - 1} \\\vdots & \quad & \quad & \quad \\1 & z_{M} & \quad & z_{M}^{Q - 1}\end{bmatrix}} & (28)\end{matrix}$

[0084] If the projection of the sets of Diracs on both xandy directionsare distinct, then rank(G)=M and the values {w₁} and {z_(i)} can beobtained from the principal left or right singular values of the matrixG, as it is well known. When there previous condition is not satisfied,the rank deficiency of G requires a different algorithm. Among the knownalternatives are the MEMP algorithm (Matrix Enhancement and MatrixPencil). The method introduces the “enhanced matrixes”, both of rank Mfrom which the sets {w₁} and {z_(i)} can be obtained, or the ACMPalgorithm (algebraic coupling of Matrix Pencils) sketched below.

[0085] Let the K(P−L)×L(N−K) enhanced matrix J be defined as$\begin{matrix}{J = \begin{bmatrix}G^{({1,1})} & G^{({2,1})} & \cdots & G^{({L,1})} \\G^{({1,2})} & \quad & \quad & \quad \\\vdots & \quad & ⋰ & \quad \\G^{({1,K})} & \cdots & \quad & G^{({L,K})}\end{bmatrix}} & (29)\end{matrix}$

[0086] where the (k,l)-th block component of J is

J_(kl)=G^((l,k)) =G _(l;P−L+l,kQ−K−1+k)  (30) $\begin{matrix}{G^{({1,k})} = \begin{bmatrix}{G\left\lbrack {l,k} \right\rbrack} & \cdots & {G\left\lbrack {l,{Q - K - 1 + k}} \right\rbrack} \\{G\left\lbrack {{l + 1},k} \right\rbrack} & \quad & {G\left\lbrack {{l + 1},{Q - K - 1 + k}} \right\rbrack} \\\vdots & \quad & \quad \\{G\left\lbrack {{P - L - 1 + l},k} \right\rbrack} & \cdots & {G\left\lbrack {{P - L - 1 + l},{Q - K - 1 + k}} \right\rbrack}\end{bmatrix}} & (31)\end{matrix}$

[0087] The matrix J can be written as

J=W′AZ′  (32)

[0088] where W′ and Z′ are generalized Vandermonde matrixes

W′=(W _(P−L) ^(T) Z _(d) W _(P−L) ^(T) Z _(d) ² W _(P−L) ^(T) . . . Z_(d) ^(K−1) W _(P−L) ^(T))  (33)

Z′=(Z _(Q−K) ^(T) W _(d) Z _(Q−K) ^(T) W _(d) ² Z _(Q−K) ^(T) . . . W_(d) ^(L−1) Z _(Q−K) ^(T))  (34)

[0089] W_(d) and Z_(d) are M×M diagonal matrixes,W_(d)=diag{exp(iω₀x_(k))}, Z_(d)=diag{exp(iω₀y_(k))} and W_(P−L) andZ_(Q−K) are given by $\begin{matrix}{W_{M - L} = \begin{bmatrix}1 & 1 & \cdots & 1 \\w_{1} & w_{2} & \quad & w_{M} \\\vdots & \quad & \quad & \quad \\w_{1}^{P - L - 1} & w_{2}^{P - L - 1} & \quad & w_{M}^{P - L - 1}\end{bmatrix}} & (35) \\{Z_{N - K} = \begin{bmatrix}1 & z_{1} & \cdots & z_{1}^{Q - K - 1} \\1 & z_{2} & \quad & z_{2}^{Q - K - 1} \\\vdots & \quad & \quad & \quad \\1 & z_{M} & \quad & z_{M}^{Q - K - 1}\end{bmatrix}} & (36)\end{matrix}$

[0090] Define the top-left matrix J_(tl) obtained by omitting the lastrow and the last column of the block components of J, i.e.$\begin{matrix}{J_{tl} = {\left. \underset{\_}{X^{({l,k})}} \right| = X_{{l:{P - L - 2 + l}},{k:{Q - K - 2 + k}}}}} & (37)\end{matrix}$

[0091] where (*)| denotes the operation of deleting the last column of(*), and (*) denotes the operation of deleting the last row of (*).Define in the similar way the top-right and bottom-left matrixes J_(tr)and J_(blt). The outline of the ACMP algorithm is then:

[0092] Compute the singular value decomposition of J_(tl)

J_(tl)=USV^(H)  (38)

[0093] Find C_(tr), C_(tb) C_(bl) and C_(br) from

U ^(H)(J _(tr) −μJ _(tl))V=F(Z _(d) −μI)G=C _(tr) −μC _(tl)  (39)

U ^(H)(J _(bl) −μJ _(tl))V=F(W _(d) −λI)G=C _(tr) −λC _(tl)  (40)

[0094] Compute the eigenvalue decomposition of the matrix

C_(tl) ⁻¹(βC _(tr)+(1−β)C _(bl))=G ⁻¹ TG  (41)

[0095]  where β is a scalar introduced with the aim of avoiding multipleeigenvalues.

[0096] Apply the eigentransformation T to find Z_(d) and W_(d)

T(C _(tl) ⁻¹ C _(bl))T ⁻¹ =W _(d)  (42)

T(C _(tl) ⁻¹ C _(tr))T ⁻¹ =Z _(d)  (43)

[0097] Since the same transformation is used to diagonalize bothmatrixes, w_(i) and z_(i) correspond to the same sinusoidal component.On the other hand, a necessary condition for this property to hold, isthat all the matrixes involved in the two matrix pencils J_(tr)−μJ_(tl)and J_(bl)−λJ_(tl) can be written as the product of the same left andright matrixes, with possibly a different matrix in the middle. Asufficient condition for KLP and Q for Y′ and Z′ to have full rank is

P−L≧M K,L≧M Q−K≧M  (44)

[0098] The relation (44) implies that if we set L=K=M, the sufficientcondition for having a unique solution is P≧2M, Q≧2M. In other words weneed to know G[m,n], 0≦m≦2M−1, 0≦m≦2N−1(in our case the G[m,n] are theFourier series coefficients of g(x,y), therefore only O(M²) samples ofg(x,y) yield a unique solution for the set {(w_(i), z_(i))}). Note thatwe assumed a deterministic signal, and as long as the condition (44) issatisfied, it is possible in this case to have a perfect reconstruction.However, in the case of noisy data, the schemes that use oversamplingyield more accurate estimation.

[0099] Finally one last step is required to solve for the correspondingset of weights {a_(i)}, which can be done starting from the equation(32), i.e. A=W′⁺JZ′⁺, where (*)⁺ denotes a pseudoinverse of (*). Sinceboth W′ and Z′ have rank M, the above system will have a unique solutionfor the matrix of coefficients A, thus the signal can be perfectlyreconstructed.

[0100] An example of the above methods can be found in FIGS. 6a to 6 d:FIG. 6a represents the original signal composed of 17 2-D Diracs, FIG.6b represents a 2-D sinc smoothing kernel, FIG. 6c represents anapproximation obtained by convolving the signal of FIG. 6a with thekernel of FIG. 6b and finally FIG. 6d shows the signal reconstructedaccording to one aspect of the invention.

[0101] In a further embodiment of the present invention, thereconstruction method of the invention is applied the case of signalmodeled as “blurred” version of sets of Diracs, or more precisely as theconvolution of the signal g(x,y) with a certain point spread function(PSF). This case is of particular interest to the fields of opticalastronomy and videogrammetry, where the image formation without noisecan be modeled as the convolution of the object containing a number ofpoint sources (i.e. stars or sharp details) with the PSF, which may be aresult of the imperfections of the imaging optics, atmospheric processesetc.

[0102] In particular the invention is used to extract the direction of adetail or a number of details from a digitized photographic image. Theinvention further allows obtaining the exact 3-D position of a number ofdetails in a set of videogrammetry images.

[0103]FIG. 4a shows a scene containing a highly contrasted point-likedetail 1. The scene is viewed by two identical cameras 10, eachproducing a 2-D projection of the scene on the respective image plane15.

[0104] The above disposition has been chosen in this example for thesake of simplifying the explanation only. The invention includes as wellmethods and devices to treat the images coming from any number ofdissimilar cameras. At the same time the scene to be analyzed does notlimit itself to a single point-like source. On the contrary any scenecontaining a plurality of details can be analyzed by the methods anddevices of the invention.

[0105] Ideally the point-like detail 1 should produce a 2-D Dirac pulsein the luminosity projected on the image plane 15. In this case it wouldbe easy to determine the exact position of detail 1, by tracing theimaginary straight line 19 joining the image 18 and the centre of theobjective 11, on which the detail 1 lays. The same is repeated for thesecond camera 20, and the two lines 19 and 29 define, with their uniqueintersection, the position of detail 1 in space.

[0106] In a real case the above procedure yields only an approximateresult. The sharpness of the image 18 is in fact degraded by the optics11, and the intensity on the image plane 15 will be described by asmooth distribution of finite width. It is usually admitted that theluminosity follows a Gaussian distribution, $\begin{matrix}{{h\left( {x,y} \right)} \propto {\exp \left( {- \frac{x^{2} + y^{2}}{2\sigma^{2}}} \right)}} & (45)\end{matrix}$

[0107] defined by σ width a describing for example the resolution of theoptical system 11. Other distribution could however be chosen, accordingto the circumstances, without leaving the scope of the claimedinvention.

[0108]FIG. 4b shows the intensity distribution on the focal plane 15 andthe distribution of the samples 50 obtained by a CCD or an equivalentimage detector 12 placed on the image plane. The finite samplingoperated by the CCD further reduces the information available, andlimits the possibility of directly locate the position of the image 18.In this case, the method above can not provide an exact position of thedetail 1, but rather an uncertainty region 100.

[0109] This situation can be improved, albeit never completely solved,by increasing the quality of the optics, reducing thus o; and adopting aCCD of higher resolution. Such solution implies however cameras ofhigher cost and size, which limits the applicability of this approach.

[0110] The light intensity coming from a finite number of point-likesources is approximated by a distribution of 2-D Dirac pulses,possessing, as discussed in the previous embodiments, a finite rate ofinnovation. The effect of the optical system 11 is to convolve theoriginal function with a smoothing function h, in the example above aGaussian kernel, while the CCD 12 performs a periodical sampling of thesmoothed signal 18. We have already shown in the previous embodimentsthat this class of signals can be perfectly reconstructed, provided asufficient number of samples are acquired. The reconstruction method ofthe invention is therefore applicable and allows exact determination ofthe directions 19, 29 of the detail 1, provided that the PSF function isknown.

[0111] In the general case, given the analytic expression for the PSF, anew sampling kernel for the invention should be h′(x,y)*φ_(g)(x,y). Thepurpose of introducing the first term h′(x,y)is to deconvolve theblurred signal prior to sampling it with the sinc sampling kernel or theGaussian sampling kernel, or another suitable kernel. Under thenoise-free assumption, this step is sufficient to obtain the set ofsamples which can be expressed as a linear combination of exponentials,and from which the location and weights of the Diracs can be uniquelyfound. However, that is true only if a is relatively small, so that theterm h′(x,y) does not blow up, otherwise the system becomesill-conditioned.

[0112] According to the circumstances, the reconstruction can beperformed directly on the 2-D data or after the application of the Radontransform, in order to reduce to one-dimensional subproblems, as in theprevious embodiments.

[0113] If more cameras are available, as in videogrammetry applications,the exact position of an object in space can be inferred by searchingthe intersection of the direction lines 19, 29 thus reconstructed, as inFIG. 4a.

[0114] The invention is not limited to the reconstruction of point-likedetails, but also to reconstruction of generic shapes in front of acontrasting background, the shapes being in general approximable with apiecewise linear or piecewise polynomial boundary. The reconstructionproceeds, in this case, by applying the appropriate differentiatedkernels ρ^((R+1)) and φ⁽²⁾ introduced in the previous embodiments.

[0115] In another embodiment of the present invention, the classicvideogrammetry problem of finding the location and/or the orientation ofone or more cameras from images of known reference points is solved. Theinvention may be employed to determine the position of a landing plane,or of another vehicle from the image of known landscape points taken bycameras on the vehicle. Also the invention may be applied to find outthe orientation of a space craft from an image of a star field taken byone or more onboard cameras. The invention may as well be applied to theproblem of determining the course of a watercraft from the images oflandscape point and/or known navigational aids.

[0116] In the general case we have C different cameras, and Ndistinguishable points visualized on each image, the position of Mpoints being known.

[0117] As it is generally known, the position and orientations in spaceof each camera are described by 6 unknowns, whereas each pointintroduces 3 unknowns. One has therefore 3·(N−4)+6·C unknowns. Theimages provide 2NC point image coordinates, one can therefore build 2NCequation where the point image coordinates figure as given values. Ifthe condition

2NC>3(N−M)+6C  (46)

[0118] is satisfied, the system is determined or overdetermined and thecamera locations and directions in space can be reconstructed.

[0119] In the case in which the locations of the N points are unknown,the problem may be solved in a relative reference system, wherein one ofthe points is placed in the origin.

[0120]FIG. 5 shows a 2-D example where one camera 120, havingcoordinates (X_(c),Z_(c)) takes one image 300 of two points 61 and 62whose known coordinates are (X_(p1),Z_(p1)) and (X_(p2),Z_(p2)).

[0121] If we determine now the coordinates v₁₁ and v₁₂ of the points 301and 302 corresponding, on the image 300, to the points 61 and 62, wehave: $\begin{matrix}\left\{ \begin{matrix}{v_{11} = {f\frac{\left( {X_{c} - X_{p1}} \right) + {{tg}\quad {\psi \left( {z_{c} - z_{p1}} \right)}}}{{{tg}\quad {\psi \left( {X_{c} - X_{p1}} \right)}} - \left( {z_{c} - z_{p1}} \right)}}} \\{v_{21} = {f\frac{\left( {X_{c} - X_{p2}} \right) + {{tg}\quad {\psi \left( {z_{c} - z_{p2}} \right)}}}{{{tg}\quad {\psi \left( {X_{c} - X_{p2}} \right)}} - \left( {z_{c} - z_{p2}} \right)}}}\end{matrix} \right. & (47)\end{matrix}$

[0122] where f represents the known focal distance.

[0123] In the particular case in which the orientation ψ of the camera120 is known and fixed, the system above becomes linear and admits oneunique solution.

[0124] The geometrical formulation of the above problems ofvideogrammetry is of course known. Actual videogrammetry applicationssuffer however from the inherent sampling limitations of real imagecapturing devices, like CCD, and from the imperfections of the opticsystem, that necessarily introduce image blurring to some degree.

[0125] The application of the sampling and reconstruction methods of thepresent invention to the above problem, provides a “superresolution”image, and allows perfect reconstruction of the position of therequested details, much beyond the limits that would be imposed by CCDand optics limitations.

[0126] By applying the reconstruction method of the invention, the givenvalues of the above problem can be determined with precision, providedthat the PSF of the cameras are known. In this way the reconstructioncan be done exactly, despite the finite resolution of the optics and ofthe CCD.

[0127] The choice of the particular geometry of the FIG. 2 was done onlyin order to exemplify the invention, and does not constitute a limitingfeature of the invention.

[0128] In the two embodiments above, reference was made to the noiselesscase and, under this hypothesis, exact reconstruction was obtained froma minimal number of samples. If however the signal has a noise componentsuperimposed, an adequate approximation of the true signal will beobtained by oversampling.

[0129] The sampling and reconstruction methods can be performed byhardware and devices in the form of dedicated circuits or systems and/orby a computer including memories in which software or firmware forcarrying out the above described methods are stored and/or executed. Theinvention relates also to a circuit, to a system and to a computerprogram for processing sets of values obtained by sampling a signal withthe above described methods.

1. Method for sampling a first multidimensional signal (g(x,y)) having afinite rate of innovation (ρ), said method comprising convoluting saidfirst signal (g(x,y)) with a sampling kernel (φ) and using a regularsampling rate (1/T_(p)), said sampling kernel and said sampling ratebeing chosen such that the sampled signal (p_(n)(θ₀)) is a completerepresentation of said first signal (g(x,y)), allowing a perfectreconstruction of said first signal, characterized in that said samplingrate (1/T_(p)) is lower than the frequency given by the Shannon theorem,but greater than or equal to the rate of innovation (p) of said firstsignal (g(x,y)).
 2. Sampling method according to the preceding claim,further comprising a preliminary step of deriving a secondone-dimensional signal (R(p,θ₀)) from said first multidimensional signal(g(x,y)).
 3. Sampling method according to claim 2, wherein saidpreliminary step of deriving said second one-dimensional signal(R(p,θ₀)) from a multidimensional signal (g(x,y)) is performed byprojecting said multidimensional signal (g(x,y)) onto a one-dimensionalline.
 4. Sampling method according to claim 20, wherein said preliminarystep of deriving said second one-dimensional signal (R(p,θ₀)) from amultidimensional signal (g(x,y)) is performed by smoothing andprojecting said multidimensional signal (g(x,y)) onto a one-dimensionalline.
 5. Sampling method according to one of the claims 3 or 4, whereinsaid step of deriving said first signal (x(t), f(p)) from amultidimensional signal (g(x,y)) is performed by applying the Radontransform to the multidimensional signal (g(x,y)) or to a smoothedversion of the multidimensional signal (g(x,y)).
 6. Sampling methodaccording to one of the claims 1 to 5, wherein said multidimensionalsignal (g(x,y)) comprises a number N of point-like features comprisingthe step of projecting said multidimensional signal (g(x,y)) onto atleast N+1 one-dimensional lines, to obtain at least N+1 one-dimensionalsignals (R(p,θ_(i))).
 7. Sampling method according to claim 6, furthercomprising the steps of: determining indicia of the said number N ofpoint-like features on each of the said N+1 one-dimensional signals(R(p,θ_(i))); determining a projecting line for each of said indicia;Inferring the position of said point-like features of saidmultidimensional signal (g(x,y)) from the positions where N+1 projectinglinecross.
 8. Sampling method according to one of the claims 1 to 7,wherein said sampling kernel (φ) is a sinc function.
 9. Sampling methodaccording to one of the claims 1 to 7, wherein said sampling kernel (φ)is a Gaussian function.
 10. Sampling method according to one of theclaims 1 to 7, wherein said sampling kernel (φ) is a spline function.11. Sampling method according to one of the claims 1 to 7, wherein saidsampling kernel (φ) is a box function.
 12. Sampling method according toone of the claims 1 to 7, wherein said sampling kernel (φ) is a hatfunction.
 13. Sampling method according to one of the claims 1 to 12,wherein said first signal (g(x,y)) is a periodic stream of weightedDirac pulses.
 14. Sampling method according to one of the claims 1 to12, wherein said multidimensional signal (g(x,y)) is a bi-level polygon.15. Sampling method according to claim 15 wherein said sampling kernel((p) is a second derivative sinc function.
 16. Sampling method accordingto one of the claims 1 to 12, wherein said multidimensional signal(g(x,y)) is a bi-level polygon with piecewise polynomial boundary with Mpieces of maximum degree R.
 17. Sampling method according to one ofclaims 1 to 12, wherein said multidimensional signal (g(x,y)) is abi-level signal with a piecewise polynomial boundary having M pieces.18. Sampling method according to claim 16 or 17, wherein said samplingkernel (φ) is a (R+1)th derivative sinc function.
 19. Method forfaithfully reconstructing a first multi-dimensional signal (g(x,y)) froma set of samples (G[m,n]), wherein the class of said signal toreconstruct (g(x,y)) is known, wherein the bandwidth of said firstmulti-dimensional signal (g(x,y)) is higher than 1T, T being thesampling interval, wherein the rate of innovation (ρ) of said faithfulreconstructed signal is finite, characterized in that said methodcomprises the step of solving a structured linear system depending onsaid known class of signal.
 20. Method according to claim 19, whereinsaid reconstruction method includes the application of singular valuedecomposition methods.
 21. Method according to claim 19, wherein saidreconstruction method includes the following steps: finding 2K spectralvalues of said first signal (g(x,y)), using an annihilating filtermethod for finding said first signal (g(x,y)) from said spectral values.22. Method according to claim 19, wherein said first signal (g(x,y)) isa finite stream of weighted Dirac pulses, said reconstruction methodincluding following steps: finding the roots of an annihilating filterto find the position of said pulses, solving a linear system to find theweights of said pulses.
 23. Method according to any of the precedingclaims, wherein said first signal (g(x,y)) has a noise componentsuperimposed, and said sampling interval (T,T_(p)) is adapted to obtainan approximation of said first signal (g(x,y)).
 24. Method according toone of the claims 1-23, wherein said first multidimensional signal(g(x,y)) is an image of a scene containing at least one feature (1, 61,62), and the exact position of said at least a feature in said image isobtained from said sampled signal (p_(n)(θ₀))
 25. Method according toclaim 24 further comprising the steps of: taking at least two images ofsaid at least one feature (1), from at least two distinct knownlocations; determining the position of said at least one feature (1) insaid scene from said exact positions in said at least two images. 26.Method according to one of the claims 1-23, wherein said firstmultidimensional signal (g(x,y)) is at least one image of a scenecontaining at least one feature (61, 62) recorded with at least oneimage recording means (120) further comprising the steps of: obtainingthe exact position of said at least one feature (61, 62) in said atleast one image; determining the position of said at least one imagerecording means (120) from said exact positions in said at least oneimage.
 27. Method according to one of the claims 1-23, wherein saidfirst multidimensional signal (g(x,y)) is at least one image of a scenecontaining at (east one feature (61, 62) recorded with at least oneimage recording means (120) further comprising the steps of: obtainingthe exact position of said at least one feature (61, 62) in said atleast one image; determining the orientation of said at least one imagerecording means (120) from said exact positions in said at least oneimage.
 28. Device comprising means operatively arranged to perform themethod of one of the preceding claims.
 29. Computer program productdirectly loadable into the internal memory of a digital processingsystem and comprising software code portions for performing the methodsof one of the preceding claims when said product is run by said digitalprocessing system.